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ABSTRACT 

It is a well-known fact that mathematical functions that are timelimited (or spacelimited) cannot be simul- 
taneously bandlimited (in frequency). Yet the finite precision of measurement and computation unavoidably 
bandlimits our observation and modeling scientific data, and we often only have access to, or are only interested 
in, a study area that is temporally or spatially bounded. In the geosciences we may be interested in spectrally 
modeling a time series defined only on a certain interval, or we may want to characterize a specific geographical 
area observed using an effectively bandlimited measurement device. It is clear that analyzing and representing 
scientific data of this kind will be facilitated if a basis of functions can be found that are "spatiospectrally" 
concentrated, i.e. "localized" in both domains at the same time. Here, we give a theoretical overview of one par- 
ticular approach to this "concentration" problem, as originally proposed for time series by Slepian and coworkers, 
in the 1960s. We show how this framework leads to practical algorithms and statistically performant methods 
for the analysis of signals and their power spectra in one and two dimensions, and on the surface of a sphere. 
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1. INTRODUCTION 

It is well appreciated that functions cannot have finite support in the temporal (or spatial) and spectral domain 
at the same time. 1 Finding and representing signals that are optimally concentrated in both is a fundamental 
problem in information theory which was solved in the early 1960s by Slepian, Landau and Pollak. 2-4 The 
extensions and generalizations of this problem 5-8 have strong connections with the burgeoning field of wavelet 
analysis. In this contribution, however, we shall not talk about wavelets, the scaled translates of a "mother" 
with vanishing moments, the tool for multi-resolution analysis. 9-11 Rather, we devote our attention entirely to 
what we shall collectively refer to as "Slepian functions" , in multiple Cartesian dimensions and on the sphere. 

These we understand to be orthogonal families of functions that are all defined on a common, e.g. geograph- 
ical, domain, where they are either optimally concentrated or within which they are exactly limited, and which 
at the same time are exactly confined within a certain bandwidth, or maximally concentrated therein. The 
measure of concentration is invariably a quadratic energy ratio, which, though only one choice out of many 12-14 
is perfectly suited to the nature of the problems we are attempting to address. These are, for example: How do 
we make estimates of signals that are noisily and incompletely observed? How do we analyze the properties of 
such signals efficiently, and how can we represent them economically? How do we estimate the power spectrum of 
noisy and incomplete data? What are the particular constraints imposed by dealing with potential-field signals 
(gravity, magnetism, etc) and how is the altitude of the observation point, e.g. from a satellite in orbit, taken 
into account? What are the statistical properties of the resulting signal and power spectral estimates? 

These and other questions have been studied extensively in one dimension, that is, for time series, but 
remarkably little work had been done in the Cartesian plane or on the surface of the sphere. For the geosciences, 
the latter two domains of application are nevertheless vital for the obvious reasons that they deal with information 
(measurement and modeling) that is geographically distributed on (a portion of) a planetary surface. In our own 
recent series of papers 15-21 we have dealt extensively with Slepian's problem in spherical geometry; asymptotic 
reductions to the plane 16,22 then generalize Slepian's early treatment of the multidimensional Cartesian case. 23 

In this chapter we provide a framework for the analysis and representation of geoscientific data by means of 
Slepian functions defined for time series, on two-dimensional Cartesian, and spherical domains. We emphasize 
the common ground underlying the construction of all Slepian functions, discuss practical algorithms, and review 
the major findings of our own recent work on signal 15 ' 17 and power spectral estimation theory on the sphere. 19, 20 



2. THEORY OF SLEPIAN FUNCTIONS 

In this section we review the theory of Slepian functions in one dimension, in the Cartesian plane, and on the 
surface of the unit sphere. The one-dimensional theory is quite well known and perhaps most accessibly presented 
in the textbook by Percival and Walden. 24 It is briefly reformulated here for consistency and to establish some 
notation. The two-dimensional planar case formed the subject of a lesser-known of Slcpian's papers 23 and is 
reviewed here also. We are not discussing alternatives by which two-dimensional Slepian functions are constructed 
by forming the outer product of pairs of one-dimensional functions. While this approach has produced some 
useful results, 25,26 it does not solve the concentration problem sensu stricto. The spherical case was treated 
in most detail, and for the first time, by ourselves elsewhere, 15-17 though two very important early studies 
by Slepian, Grunbaum, and others, laid much of the foundation for the analytical treatment of the spherical 
concentration problem for cases with special symmetries. 27,28 Finally, we recast the theory in the context of 
reproducing- kernel Hilbert spaces, through which the reader may appreciate some of the connections with radial 
basis functions, splines, and wavelet analysis, which are commonly formulated in such a framework. 29 

2.1. Spatiospectral concentration for time series 
General theory in one dimension 

We use t to denote time or one-dimensional space and u for angular frequency, and adopt a normalization 
convention 11 in which a real-valued time-domain signal f(t) and its Fourier transform F(oj) are related by 

/oo />oo 
FHe^dw, F(oj)= / f(t)e^ l dt. (1) 
-oo J — oo 

The problem of finding the strictly bandlimited signal 

/w 
G{u)e^ dej, (2) 
-w 

that is maximally (though by virtue of the Paley- Wiener theorem 9, 11 never completely) concentrated into a time 
interval \t\ < T was first considered by Slepian, Landau and Pollak. 2 ' 3 The optimally concentrated signal is 
taken to be the one with the least energy outside of the interval: 

£ T 9 2 (t)dt 

A = — -Too = maximum. (3) 

/ 9 2 (t)dt 

J — oo 

Bandlimited functions g{t) satisfying the variational problem (3) have spectra G(uj) that satisfy the frequency- 
domain convolutional integral eigenvalue equation 



/ 



w 

D(u, w') G{uj') duj' = XG(uj), \u\ < W, (4a) 



w 



tt(lu — U)') 

The corresponding time- or spatial-domain formulation is 

i-T 



D(w,w') = 7^ jr 1 - ( 4b ) 



J ^D(t,t')g(t')dt' = Xg(t), 46R, 



(5a) 



(5b) 



The "prolate spheroidal cigentapers" gi(t), g2(t), ■ ■ ■ that solve eq. (5) form a doubly orthogonal set. When they 
are chosen to be orthonormal over infinite time \t\ < oo they are also orthogonal over the finite interval \t\ < T: 



/oo pT 
g a gp dt = 5 a p, / g a gp dt = \ a S af3 . (6) 
-oo J-T 



rT 
-T 

A change of variables and a scaling of the eigenfunctions transforms eq. (4) into the dimcnsionless eigenproblcm 



/ 



l 

D(x,x')ip(x')dx' = \ip{x), (7a) 



smTW(x-x') 

D(x, x 1 ) = . 7b 

7T(X — X ) 

Eq. (7) shows that the eigenvalues Ai > A 2 > . . . and suitably scaled eigenfunctions ipi(x), ^2 0*0, ■ ■ ■ depend only 
upon the time-bandwidth product TW. The sum of the concentration eigenvalues A relates to this product by 



oo „i 
N 1D = A « = / D(x,x) 

a=l 



dx = (2T)(2W) = 2TW 

27T 7T W 



The shape of the eigenvalue spectrum of eq. (7) has a characteristic step shape, showing significant (A s=s 1) and 
insignificant (A w 0) eigenvalues separated by a narrow transition band. 30 ' 31 Thus, this "Shannon number" is 
a good estimate of the number of significant eigenvalues, or, roughly speaking, iV 1D is the number of signals 
f(t) that can be simultaneously well concentrated into a finite time interval \t\ < T and a finite frequency 
interval |w| < W. In other words, 4 iV 1D is the approximate dimension of the space of signals that is "essentially" 
timelimited to T and bandlimited to W, and using the orthogonal set gi, c/2, ■ • ■ , gN 1D as its basis is parsimonious. 

Sturm-Liouville character and tridiagonal matrix formulation 

The integral operator acting upon ip in eq. (7) commutes with a differential operator that arises in expressing 
the three-dimensional scalar wave equation in prolate spheroidal coordinates, 1, 2 which makes it possible to find 
the scaled eigenfunctions tjj by solving the Sturm-Liouville equation 



d_ 

dx 



2^dip 



+ 



x- 



(iV 1D ) 2 7T 2 2 



1> = 0, (9) 



where x 7^ A is the associated eigenvalue. The eigenfunctions ij){x) of eq. (9) can be found at discrete values of x 
by diagonalization of a simple symmetric tridiagonal matrix 24 ' 28 ' 32 with elements 

{[N - 1 - 2x]/2) 2 cos(27rVT) for x = 0, • • • , N - 1, 

x(N-x)/2 for ar = l,...,iV-l. (10) 

The matching eigenvalues A can then be obtained directly from eq. (7). The discovery of the Sturm-Liouville 
formulation of the concentration problem posed in eq. (3) proved to be a major impetus for the widespread 
adoption and practical applications of the "Slepian" basis in signal identification, spectral analysis and numerical 
analysis. Compared to the sequence of eigenvalues A, the spectrum of the eigenvalues \ is extremely regular and 
thus the solution of eq. (9) is without any problem amenable to finite-precision numerical computation. 24 

2.2. Spatiospectral concentration in the Cartesian plane 
General theory in two dimensions 

A square- integrable function /(x) defined in the plane has the two-dimensional Fourier representation 



/oo poo 
F(k)e 4kx dk, F(k)= / .f(x) e - lkx dx, 
-OO J — OO 



(11) 



We use g(x) to denote a function that is bandlimited to AC, an arbitrary subregion of spectral space, 

ff(x) = (2tt)- 2 / G(k)e lkx dk. (12) 



Following Slepian, 23 we seek to concentrate the power of g(x) into a finite spatial region 1Z e R 2 , of area A: 

/ ff 2 (x)dx 

A = -753 = maximum. (13) 

/ 5 2 (x)dx 

J — oo 

Bandlimited functions g(x) that maximize the Rayleigh quotient (13) solve the Frcdholm integral equation 

L>(k, k') G(k') dk' = AG(k), k e AC, (14a) 



Ik 

D(k,k') - (2vr)- 2 / e 4 ( k '- k ) x dx. 
The corresponding problem in the spatial domain is 



(14b) 



D(x,x')g(x')dx' = Xg(x), xeR 2 , (15a) 



£>(x,x') = (2^)- 2 / e lk (x - x,) dk. 
Jk 



(15b) 



The bandlimited spatial-domain eigenfunctions gi(x), <y 2 (x), . . . and eigenvalues Ai > A 2 > . . . that solve eq. (15) 
may be chosen to be orthonormal over the whole plane ||x|| < oo in which case they are also orthogonal over 1Z: 



/ g a g(i dx = S a p, / g a gp dx = X a S a0 . 



(16) 



Concentration to the disk-shaped spectral band AC = {k : ||k|| < K} allows us to rewrite eq. (15) after a change 
of variables and a scaling of the eigenfunctions as 

/ £(£,OV>(£V£' = AV(£), (17a) 
Jn, 

2tt ' 

where the region TZ* is scaled to area An and J\ is the first-order Bessel function of the first kind. Eq. (17) shows 
that, also in the two-dimensional case, the eigenvalues Ai, A 2 , . . . and the scaled eigenfunctions ip2(£), ■ ■ ■ 

depend only on the combination of the circular bandwidth K and the spatial concentration area A, where the 
quantity K 2 A/(4ir) now plays the role of the time-bandwidth product TW in the one-dimensional case. The 
sum of the concentration eigenvalues A defines the two-dimensional Shannon number iV 2D as 

» 2D = j> = jU«.e« - - a- A m 

Just as N 1B in eq. (8), N 2B is the product of the spectral and spatial areas of concentration multiplied by the 
"Nyquist density". 5 ' 9 And, similarly, it is the effective dimension of the space of "essentially" space- and band- 
limited functions in which the set of two-dimensional functions g ll g 2l ■ ■ ■ , gN 2D may act as a sparse orthogonal 
basis. 

An example of Slepian functions on a geographical domain in the Cartesian plane can be found in Figure 1. 
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Figure 1. Bandlimited eigenfunctions g\, g2, . . . , Qa that are optimally concentrated within the Columbia Plateau, a 
physiographic region in the United States centered on 116.02°W 43.56°N (near Boise City, Idaho) of area A ~ 145 x 
10 3 km 2 . The concentration factors Ai, A2, . . . , A4 are indicated; the Shannon number iV 2D = 10. The top row shows a 
rendition of the eigenfunctions in space on a grid with 5 km resolution in both directions, with the convention that positive 
values are blue and negative values red, though the sign of the functions is arbitrary. The spatial concentration region is 
outlined in black. The bottom row shows the squared Fourier coefficients |G a (k)| 2 as calculated from the functions g a (x) 
shown, on a wavenumber scale that is expressed as a fraction of the Nyquist wavenumber. The spectral limitation region 
is shown by the black circle at wavenumber K = 0.0295 rad/km. All areas for which the absolute value of the functions 
plotted is less than one hundredth of the maximum value attained over the domain are left white. The calculations were 
performed by the Nystrom method using Gauss-Legendre integration of eq. (17) in the two-dimensional spatial domain. 22 

Sturm-Liouville character and tridiagonal matrix formulation 

If in addition to the circular spectral limitation, space is also circularly limited, in other words, if the spatial 
region of concentration or limitation 1Z is a circle of radius i?, then a polar coordinate, x = (r, 9), representation 



y/2g(r) cos m9 if m < 0, 
9{r, °) = { g{ r ) if m = 0, 

\[2 g(r) smm9 if m > 0, 



(19) 



may be used to decompose eq. (17) into a series of non-degenerate fixed-order eigenvalue problems, after scaling, 



47V / J m (2VN™p£)J m (2VN™pZ')pdp. 



The solutions to eq. (20) also solve a Sturm-Liouville equation on < £ < 1. In terms of </?(£) = \/f 



dip 



(20a) 
(20b) 

(21) 



for some A. When m = ±1/2 eq. (21) reduces to eq. (9). By extension to £ > 1 the fixed-order "generalized 
prolate spheroidal functions" <£i(£), <^2(0> • • ■ can be determined from the rapidly converging infinite series 



^ ^ (i + my. V c £ 



(22) 



where the fixed-m expansion coefficients di are determined by recursion 23 or by diagonalization of a non- 
symmetric tridiagonal matrix 33,34 with elements given by 



71+1/ 



c 2 (m + l + l) 2 



(2l + m+ l)(2Z + m + 2)' 



T H = 2? + m 



Th+i — - 



2 

c 2 (Z + l) 2 



2Z + TO + - + 



(2Z + m)(2Z + m + 2) 



(2/ + m + 2)(2Z + m + 3)' 



(23) 



where the parameter Z ranges from to some large value that ensures convergence. The desired concentration 
eigenvalues A can subsequently be obtained by direct integration of eq. (17), or, alternatively, from 



A = 2 7 V7V 2D , with 7 



2 ™+i(m + l)! ^ < y 



(24) 



Numerous numerical methods exist to use eqs (14)-(15) and (21) in solving the concentration problem (13) 
for a variety of two-dimensional Cartesian domains. After a long hiatus since the work of Slepian, 23 the two- 
dimensional problem has recently been the focus of renewed attention in the applied mathematics community, 33 ' 34 
and applications to the geosciences are to follow 22 . 

An example of Slepian functions on a disk-shaped region in the Cartesian plane can be found in Figure 2. 

2.3. Spatiospectral concentration on the surface of a sphere 
General theory in "three" dimensions 

We denote the colatitude of a geographical point f on the unit sphere surface = {f : ||?|| = 1} by < 9 < tt and 
the longitude by < (f> < 2ir. We use R to denote a region of f2, of area A, within which we seek to concentrate 
a bandlimitcd function of position f. We express a square- integrable function /(f) on the surface of the unit 
sphere as 



/( r ) = Y flm,Yl m {r), flm 

1=0 m=-l 

using orthonormalized real surface spherical harmonics 35,36 

V2X lm (9)cos 



Y lm (r) = Y lm (e, 



Xl m (6) 
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V2X lm (6)s\nmq 

(-1) 
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erf 
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(l + m)\_ 



f(r)Y lm (r)dn, 



if -I < m < 0, 
if m = 0, 
if < m < I, 

1/2 



(25) 



P /m (cosfl), 



2 l U 
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Figure 2. Bandlimited eigenfunctions g a (r,8) that are optimally concentrated within a Cartesian disk of radius R=l. 
The dashed circle denotes the region boundary. The Shannon number iV 2D = 42. The eigenvalues X a have been sorted 
to a global ranking with the best concentrated eigenfunction plotted at the top left and the 30th best in the lower right. 
Blue is positive and red is negative and the color axis is symmetric, but the sign is arbitrary; regions in which the absolute 
value is less than one hundredth of the maximum value on the domain are left white. The calculations were performed 
by Gauss-Legendre integration in the two-dimensional spatial domain, which sometimes leads to slight differences in the 
last two digits of what should be identical eigenvalues for each pair of non-circularly-symmetric eigenfunctions. 



The quantity < I < oo is the angular degree of the spherical harmonic, and — I < m < I is its angular order. 
The function Pi m {p) defined in (28) is the associated Legendre function of integer degree I and order m. Our 
choice of the constants in eqs (26)-(27) orthonormalizes the harmonics on the unit sphere: 



/ 



Yi m Yi> m > dCl = 5w5 mm >, (29) 



in 

and leads to the addition theorem in terms of the Legendre functions Pi (fi) = Pi (fi) as 

l 



J2 YUmrn(r')=(^)p l (r-r'). 

m= — l ^ ' 



(30) 



To maximize the spatial concentration of a bandlimited function 

l i 

•9( f ) = E E 9imYi m (i) (31) 
within a region R, we maximize the energy ratio 

/ 9 2 (r)dn 

\ - Jr . = maximum. (32) 



g 2 (r) dQ 

n 



Maximizing equation (32) leads to the positive-definite spectral-domain eigenvalue equation 

L I' 

E E D lm,l'm'9l'm' = ^9lm, < I < L, (33a) 

Z'=0m'=-Z' 



Am.i'm' = / Yi m Y Vm > dfl, (33b) 

JR 

and we may equally well rewrite eq. (33) as a spatial-domain eigenvalue equation: 

D(r,r')g{r')dQ' = \g(r), fefi, (34a) 



^ f, ) = E(^) p '( f - f ')' ( 34b ) 

where Pi is the Legendre function of degree I. Eq. (34) is a homogeneous Frcdholm integral equation of the second 
kind, with a finite-rank, symmetric, Hermitian kernel. We choose the spectral eigenfunctions of the operator in 
eq. (33b), whose elements are gi ma , a = l, . . . , (L + l) 2 , to satisfy the orthonormality relations 

L L L 

E^imaS/m/S = Sa/3, E 9 lm a E ^ImA' m' 9V "m' = A Q (5 a /3. (35) 

lm Ira I'm' 

The finite set of bandlimited spatial eigensolutions ffi(f ), 32(f), ••• , 5(l+i) 2 (?) can be made orthonormal over the 
whole sphere SI and orthogonal over the region R: 

/ g a gpdfl^ Sap, / g a gpdQ = X a 5 a p. (36) 

In the limit of a small area A — > and a large bandwidth L — > 00 and after a change of variables, a scaled version 
of eq. (34) will be given by 

/ D(£,€)M)d£t. = \1>{t), (37a) 

JR. 



where the scaled region i?» now has area A-k and J\ again is the first-order Bessel function of the first kind. 
As in the one- and two-dimensional case, the asymptotic, or "flat-Earth" eigenvalues Ai > A 2 > . . . and scaled 
eigenfunctions ipi(£), ip2(£), ■ ■ ■ depend upon the maximal degree L and the area A only through what is once 
again a space-bandwidth product, the "spherical Shannon number", this time given by 



iV 3D = A «=E E D lm,lm= f D(T,T)dQ 

a=l 1=0 m=-l J R 

= f ^,$)dQ t = (L+l) 2 ^-. (38) 

JR, 47r 

Irrespectively of the particular region of concentration that they were designed for, the complete set of band- 
limited spatial Slepian eigenfunctions <7i,<72> • • ■ ; .9(l+i) 2 j is a basis for bandlimited scalar processes anywhere on 
the surface of the unit sphere. 16 ' 17 This follows directly from the fact that the spectral localization kernel (33b) 
is real, symmetric, and positive definite: its eigenvectors gn m , giim, ■ ■ ■ , 9(L+i) 2 im form an orthogonal set as we 
have seen. Thus the Slepian basis functions <?«(?), a = 1, . . . , (L + l) 2 given by eq. (31) simply transform the 
same-sized limited set of spherical harmonics Yi m (r), < I < L, —l<m<l that are a basis for the same space 
of bandlimited spherical functions with no power above the bandwidth L. The effect of this transformation is 
to order the resulting basis set such that the energy of the first iV 3D functions, gi(r), . . . ,g N 3n(r), with eigen- 
values Awl, are concentrated in the region R, whereas the remaining eigenfunctions, gjv 3D +i(f); • ■ ■ 7 9{l+i) 2 (r ), 
are concentrated in the complimentary region O — R. As in the one- and two-dimensional case, therefore, the 
reduced set of basis functions gi,g2, ■ ■ ■ , <?./v 3D can be regarded as a sparse, global, basis suitable to approximate 
bandlimited processes that are primarily localized to the region R. The dimensionality reduction is dependent 
on the fractional area of the region of interest. In other words, the full dimension of the space (L + l) 2 can be 
"sparsified" to an effective dimension of N 3B = (L + l) 2 ^4/(47r) when the signal of interest resides in a particular 
geographic region. 

An example of Slepian functions on a geographical domain on the surface of the sphere is found in Figure 3. 
Sturm-Liouville character and tridiagonal matrix formulation 

In the special but important case in which the region of concentration is a circularly symmetric cap of colatitudinal 
radius <d, centered on the North Pole, the colatitudinal parts g{9) of the separable functions 

y/2 g (6) cos m(f> if — L < m < 0, 
<l(B.n) = { g (9) if m = 0, (39) 

V2g(9) sin m<p if < m < L, 



which solve eq. (34), or, indeed, the fixed-order versions 
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,e 

/ D(6, 6')g{6') sin 6' d6' = \g{6), < 9 < tt, (40a) 



D{9, 9') = 2nJ2 X lm (6)X lm (6'), (40b) 



are identical to those of a Sturm-Liouville equation for the g{9). In terms of ^ = cos 9 
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Figure 3. Bandlimited L = 60 eigenfunctions gi, g%, . . . , gi2 that are optimally concentrated within Antarctica. The 
concentration factors Ai, A2, • • • , A12 are indicated; the rounded Shannon number is iV 3D = 102. The order of concentration 
is left to right, top to bottom. Positive values are blue and negative values are red; the sign of an eigenfunction is arbitrary. 
Regions in which the absolute value is less than one hundredth of the maximum value on the sphere are left white. We 
integrated eq. (33b) over a splined high-resolution representation of the boundary, using Gauss-Legendre quadrature over 
the colatitudes, and analytically in the longitudinal dimension. 18 



with x 7^ A. This equation can be solved in the spectral domain by diagonalization of a simple tridiagonal matrix 
with a very well-behaved spectrum. 16, 18 This matrix, whose eigenfunctions correspond to the gi m of eq. (31) at 
constant m is given by 



T u = -1(1 + 1) cosQ, 
T ll+ i = [1(1 + 2) - L(L + 2)] , 



(2/ + 1)(2Z + 3) ' 



(42) 



Moreover, when the region of concentration is a pair of axisymmetric polar caps of common colatitudinal radius 
centered on the North and South Pole, the g(9) can be obtained by solving the Sturm-Liouville equation 



d 



( M 2-cos 2 e)(l- M 2 )^ 



+ [x + L p (L p + 3)n J - 



m 2 (p 2 - cos 2 6) 
1~7 



.9 = 0, 



(43) 



where L p = L or L p = L — 1 depending whether the order m of the functions g(9,<f>) in eq. (39) is odd or 
even, and whether the bandwidth L itself is odd or even. In their spectral form the coefficients of the optimally 
concentrated antipodal polar-cap eigenfunctions only require the numerical diagonalization of a tridiagonal matrix 
with analytically prescribed elements and a spectrum of eigenvalues that is guaranteed to be simple, 17, 18 namely 

1% = -l(l + l)cos 2 e + ^[(l + l) 2 -m 2 ] 

+ [(l-2)(l + l)-L p (L p + ^ 1 2 3 - 2 -^ + D 



3 3 (2Z + 3)(2Z- 1) 



[1(1 + 3) - L P (L P + 3)] / [(/ + 2) 2 - m 2 ] [{I + 1) 



m 2 



T "+ 2 21 + 3 V (2/ + 5)(2/ + 1) ' (44) 

The concentration values A, in turn, can be determined from the defining equations (33) or (34). The spectra 
of the eigenvalues % of eqs (42) and (44) display roughly equant spacing, without the numerically troublesome 
plateaus of nearly equal values that characterizes the eigenvalues A. Thus, for the special cases of symmetric 
single and double polar caps, the concentration problem posed in eq. (32) is not only numerically feasible also 
in circumstances where direct solution methods are bound to fail, 37 but essentially trivial in every situation. 
Numerical methods for the solution of eqs (33)-(34) on completely general domains on the surface of the sphere 
were discussed by us elsewhere. 16-18 

2.4. Mid-term summary 

It is interesting to reflect, however heuristically, on the commonality of all of the above aspects of spatiospectral 
localization, in the slightly expanded context of rcproducing-kcrncl Hilbert spaces. 9 ' 38-40 In one dimension, the 
Fourier orthonormality relation and the "reproducing" properties of the spatial delta function are given by 

/CO POO 
e^-t'Uw, / f(t')S(t,t')dt' = f(t). (45) 
-oo J —CO 

In two Cartesian dimensions the equivalent relations are 

/OO /'CO 
e* k (x-x') dK / /(x ' )(5(x , x ') dx ' = /(x ), (46) 
-co J —oo 

and on the surface of the unit sphere we have 

5(f,f') = £(^)fl(f-f'), Jj(r')S(r,r')dn' = f(r). (47) 

The integral-equation kernels (5b), (15b) and (34b) are all bandlimitcd spatial delta functions which are repro- 
ducing kernels for bandlimited functions of the types in eqs (2), (12) and (31): 



/VV />OC 
e^-^duj, / g(t')D(t,t')dt' = g(t), (48) 

-W J-oo 

D(x,x') - (2tt)- 2 I e 4k ( x - x ') dk, [ .g(x')L'(x,x , )dx' = g(x), (49) 

J K, J — oo 

D(f,f') = ^(^l)P;(f-f'), lj(r')D(r,r / )dn = g(r). (50) 

The equivalence of eq. (48) with eq. (5b) is through the Euler identity and the reproducing properties follow from 
the spectral forms of the orthogonality relations (45)-(46), which are self-evident by change-of- variables, and from 
the spectral form of eq. (47) , which is eq. (29) . Much as the delta functions of eqs (45)-(47) set up the Hilbert 
spaces of all square-integrable functions on the real line, in two-dimensional Cartesian space, and on the surface 



of the sphere, the kernels (48)— (50) induce the equivalent subspaces of bandlimited functions in their respective 
dimensions. Inasmuch as the Slepian functions are the integral eigenfunctions of these reproducing kernels in the 
sense of eqs. (5a), (15a) and (34a), they are complete bases for their hand-limited subspaces. 2,3 ' 9,10,29 Therein, 
the N 1T> , 7V 2D or A^ 3D best time- or space- concentrated members allow for sparse, approximate, expansions of 
signals that are spatially concentrated to the one-dimensional interval t G [— T, T] c R, the Cartesian region 
x G 1Z C R 2 , or the spherical surface patch f G R C O. 

As a corollary to this behavior, the infinite sets of exactly time- or spacelimited (and thus band- concentrated) 
versions of the functions g, which are the eigenfunctions of eqs (5), (15) and (34) with the domains appropri- 
ately restricted, are complete bases for square-integrable scalar functions on the intervals to which they are 
confined. 2 ' 3, 16 Expansions of such wideband signals in the small subset of their N ir> , N 2B or 7V 3D most band- 
concentrated members provide reconstructions which are constructive in the sense that they progressively capture 
all of the signal in the mean squared sense, in the limit of letting their numbers grow to infinity. This second class 
of functions can be trivially obtained, up to a multiplicative constant, from the bandlimited Slepian functions g 
by simple time- or space limitation. While Slepian, 1,2 for this reason perhaps, never gave them a name, we have 
been referring to those as h in our own investigations of Slepian functions on the sphere. 16,17,20 



3. PROBLEMS IN THE GEOSCIENCES AND BEYOND 

Taking all of the above at face-value but referring again to the literature cited thus far for proof and additional 
context, we return to considerations closer to home, namely the estimation of geophysical (or cosmological) 
signals and/or their power spectra, from noisy and incomplete observations collected at or above the surface 
of the spheres "Earth" or "planet" (or from inside the sphere "sky"). We restrict ourselves to real- valued 
scalar measurements, contaminated by additive noise for which we shall adopt idealized models. We focus 
exclusively on data acquired and solutions expressed on the unit sphere. We have considered generalizations 
to problems involving satellite data collected at an altitude and/or potential fields elsewhere. 15,17 ' 18 ' 20 Two 
different statistical problems come up in this context, namely, (i) how to find the "best" estimate of the signal 
given the data, and (ii) how to construct from the data the "best" estimate of the power spectral density of the 
signal in question. 

Thus, let there be some data distributed on the unit sphere, consisting of "signal" , s and "noise" , n, and let 
there be some region of interest R C fi, in other words, let 

dm = { s ( f ) + n ( f ) if f g i?, 

^ ' [ unknown/undesired if f G R — Q 

We assume that the signal of interest can be expressed by way of spherical harmonic expansion as in eq. (25), 
and that it is, itself, a realization of a zero-mean, Gaussian, isotropic, random process, namely 

oo / 

s ( f )=X] ^2 s lm Y lm (r), (sim) = and {si m s Vm >) = Si S w d mm >. (52) 

1=0 m=-l 

For illustration we furthermore assume that the noise is a zero-mean stochastic process with an isotropic power 
spectrum, i.e. (n(r)) = and (ni m ni' m i) — Ni 5u'5 mm > , and that it is statistically uncorrelated with the signal. 
We refer to power as white when Si = S or Ni = N, or, equivalently, when (n(f)n(f')) = NS(r, f). Our objectives 
are thus (i) to determine the best estimate §i m of the spherical harmonic expansion coefficients s; TO of the signal 
and (ii) to find the best estimate Si for the isotropic power spectral density Si . While in the physical world there 
can be no limit on bandwidth, practical restrictions force any and all of our estimates to be bandlimited to some 
maximum spherical harmonic degree L, thus by necessity §i m = and Si = for I > L: 

l i 

= £ *Jm*Jm(f). (53) 

1=0 m=-l 

This limitation, combined with the statements eq. (51) on the data coverage or the study region of interest, 
naturally leads us back to the concept of "spatiospectral concentration", and, as we shall see, solving either 



(51) 



problem (i) or (ii) will gain from involving the "localized" Slepian functions rather than, or in addition to, the 
"global" spherical harmonics basis. 

This leaves us to clarify what we understand by "best" in this context. While we adopt the traditional 
statistical metrics of bias, variance, and mean squared error to appraise the quality of our solutions, 41 ' 42 the 
resulting connections to sparsity will be real and immediate, owing to the Slepian functions being naturally 
instrumental in constructing efficient, consistent and/or unbiased estimates of either s; m or Si. Thus, we define 

v = (s 2 )-(s) 2 , b=(s)-s, e = s-s, and (e 2 ) = v + b 2 (54) 

for problem (i), where the lack of subscript indicates that we can study variance, bias and mean squared error of 
the estimate of the coefficients §i m but also of their spatial expansion s(r). For problem (ii) on the other hand, 
we focus on the estimate of the isotropic power spectrum at a given spherical harmonic degree I by identifying 

vi = (Sf) - (Si) 2 , bi = (Si)-Si, ei = Si-Si, and (e 2 ) = v t + b 2 . (55) 

Depending on the application, the "best" estimate could mean the unbiased one with the lowest variance, 43-47 it 
could be simply the minimum- variance estimate having some acceptable and quantifiable bias, 19 or, as we would 
usually prefer, it would be the one with the minimum mean squared error. 17 ' 20 

3.1. Problem (i): Signal estimation from noisy and incomplete spherical data 
Spherical harmonic solution 

Paraphrasing results elaborated elsewhere, 17 we write the bandlimited solution to the damped inverse problem 



(s — d) 2 dCl + rj I s 2 dfl — minimum, 

Jr. Jr. 



(56) 



where rj > is a damping parameter, by straightforward algebraic manipulation, as 

L 



Sim = ^ ^2 ( D lm,l'm' + V D lm,l'm') / dY Vm > d£l, (57) 



l'=0 m' = -l 

where Dim,Vm'^ the kernel that localizes to the region Q — R, compliments -Dz m ,Z'm' given by eq. (33b) which 
localizes to R. Given the eigenvalue spectrum of the latter, its inversion is inherently unstable, thus eq. (56) is an 
ill-conditioned inverse problem unless 77 > 0, as has been well known, e.g. in geodesy. 48 ' 49 Elsewhere 17 we have 
derived exact expressions for the optimal value of the damping parameter r\ as a function of the signal-to-noise 
ratio under certain simplifying assumptions. As can be easily shown, without damping the estimate is unbiased 
but effectively incomputable; the introduction of the damping term stabilizes the solution at the cost of added 
bias. And of course when R = ft, eq. (57) is simply the spherical harmonic transform, as in that case, eq. (33b) 
reduces to eq. (29), in other words, then £>; TO; ;' m ' = Su<5 mm <. 

Slepian basis solution 

The trial solution in the Slepian basis designed for this region of interest R, i.e. 

(L+lf 

= § c9a(r), (58) 

a=l 

would be completely equivalent to the expression in eq. (53) by virtue of the completeness of the Slepian basis 
for bandlimited functions everywhere on the sphere and the unitarity of the transform (31) from the spherical 
harmonic to the Slepian basis. The solution to the undamped (77 = 0) version of eq. (56) would then be 



= A" 1 / 

JR 



dg a dil, (59) 



which, being completely equivalent to eq. (57) for rj = 0, would be computable, and biased, only when the 
expansion in eq. (58) were to be truncated to some finite J < (L + 1) 2 to prevent the blowup of the eigenvalues A. 
Assuming for simplicity of the argument that J = N 3B , the essence of the approach is now that the solution 

JV 3D 

«(r) = ^2 s a g a (r) (60) 

a=l 

will be sparse (in achieving a bandwidth L using 7V 3D Slepian instead of (L + l) 2 spherical-harmonic expansion 
coefficients) yet good (in approximating the signal as well as possible in the mean squared sense in the region of 
interest R) and of geophysical utility (assuming we are dealing with spatially localized processes that are to be 
extracted, e.g., from global satellite measurements). 21 ' 50 

Bias and variance 

In concluding this section let us illustrate another welcome by-product of our methodology, by writing the mean 
squared error for the spherical harmonic solution (57) compared to the equivalent expression for the solution in 
the Slepian basis, eq. (59). We do this as a function of the spatial coordinate, in the Slepian basis for both, and, 
for maximum clarity of the exposition, using the contrived case when both signal and noise should be white as 
well as bandlimited (which technically is impossible). In the former case, we get 

(L + l) 2 

(e 2 (f)) = N X a [X a+V (l-X a )]- 2 gl(r) 

a=l 

(L+l) 2 

+ v 2 s J2 (i-K) 2 [K + v(i-KT 2 g 2 a (t), 

a=l 

while in the latter, we obtain 

N 3D (L + l) 2 

(e 2 (r))=Nj2K 1 9l(?) + S £ <£(f). (62) 

a=l a>JV 3D 

All (L+ l) 2 basis functions are required to express the mean squared estimation error, whether in eq. (61) or in 
eq. (62). The first term in both expressions is the variance, which depends on the measurement noise. Without 
damping or truncation the variance grows without bounds. Damping and truncation alleviate this at the expense 
of added bias, which depends on the characteristics of the signal, as given by the second term. In contrast to 
eq. (61), however, the Slepian expression (62) has disentangled the contributions due to noise/ variance and 
signal/bias by projecting them onto the sparse set of well-localized and the remaining set of poorly localized 
Slepian functions, respectively. The estimation variance is felt via the basis functions a = 1 — ► N 3I> that are well 
concentrated inside the measurement area, and the effect of the bias is relegated to those a — N 3B + 1 — > (L + l) 2 
functions that are confined to the region of missing data. 

When forming a solution to problem (i) in the Slepian basis by truncation according to eq. (60), changing the 
truncation level to values lower or higher than the Shannon number N 3B amounts to navigating the trade-off 
space between variance, bias (or "resolution"), and sparsity in a manner that is captured with great clarity by 
cq. (62). We refer the reader elsewhere 17, 18 for more details, and, in particular, for the case of potential fields 
estimated from data collected at satellite altitude. 

3.2. Problem (ii): Power spectrum estimation from noisy and incomplete spherical data 

Following 20 we will find it convenient to regard the data d(r) given in eq. (51) as having been multiplied by a 
unit-valued boxcar window function confined to the region R, 

p=0q=—p k 



(61) 



The power spectrum of the boxcar window (63) is 



2p + 



1 p 
— Tb 2 



(64) 



The spherical periodogram 

Should we decide that an acceptable estimate of the power spectral density of the available data is nothing else 
but the weighted average of its spherical harmonic expansion coefficients, we would be forming the spherical 
analogue of what Schuster 51 named the "periodogram" in the context of time series analysis, namely 



J r-U)5TT£[/,«" 1 W'>« 



(65) 



Bias of the periodogram 

Upon doing so we would discover that the expected value of such an estimator would be the biased quantity 



(Sf p ) = J2K w (S l ,+N l , 



(66) 



where, as it is known in astrophysics and cosmology 



=o 

52-54 



the periodogram "coupling" matrix 
l l' 



KW = 2TTT S S i D lm,l'm'} 

m— — l rn' — — l' 



(67) 



governs the extent to which an estimate Sf F of Si is influenced by spectral leakage from power in neighboring 
spherical harmonic degrees I' = I ± 1, 1 ± 2, . . ., all the way down to and up to oo. In the case of full data 
coverage, R = f2, or of a perfectly white spectrum, Si — S, however, the estimate would be unbiased — provided 
the noise spectrum, if known, can be subtracted beforehand. 

Variance of the periodogram 

The covariance of the periodogram estimator (65) would moreover be suffering from strong wideband coupling 
of the power spectral densities in being given by 



v sp 



2(4tt/A) 2 



(2/ + 1)(2Z' + 1) 



E E 



-i 2 



m' = — i' Lp=0 <j=0 



(68) 



Even under the commonly made assumption as should the power spectrum be slowly varying within the main 
lobe of the coupling matrix, such coupling would be nefarious. In the "locally white" case we would have 



V SP 



(2TT1WT1) & + N l){ S v +Nv) E. EJ^.-] 2 - 



m— — l m' =—V 

Only in the limit of whole-sphere data coverage will eqs (68) or (69) reduce to 



v w 



WS 



2 (St + N^Sa^ 



21 + 1 



(69) 



(70) 



which is the "planetary" or "cosmic" variance that can be understood on the basis of elementary statistical 
considerations. 55-57 The strong spectral leakage for small regions (A <g; An) is highly undesirable and makes the 
periodogram 'hopelessly obsolete', 58 or, to put it kindly, 'naive', 24 just as it is for one-dimensional time series. 



In principle it is possible — after subtraction of the noise bias — to eliminate the leakage bias in the 
periodogram estimate (65) by numerical inversion of the coupling matrix Kw . Such a "deconvolved periodogram" 
estimator is unbiased. However, its covariance depends on the inverse of the periodogram coupling matrix, which 
is only feasible when the region R covers most of the sphere, A w Aw. For any region whose area A is significantly 
smaller than Att, the periodogram coupling matrix (67) will be too ill-conditioned to be invertible. 

Thus, much like in problem (i) we are faced with bad bias and poor variance, both of which are controlled 
by the lack of localization of the spherical harmonics and their non-orthogonality over incomplete subdomains 
of the unit sphere. Both effects are described by the spatiospectral localization kernel defined in (33b), which, in 
the quadratic estimation problem (ii) appears in "squared" form in eq. (68) . Undoing the effects of the wideband 
coupling between degrees at which we seek to estimate the power spectral density by inversion of the coupling 
kernel is virtually impossible, and even if we could accomplish this to remove the estimation bias, this would 
much inflate the estimation variance. 20 

The spherical multitaper estimate 

We therefore take a page out of the one-dimensional power estimation playbook of Thomson 59 by forming 
the "eigenvalue-weighted multitaper estimate" . We could weight single-taper estimates adaptively to minimize 
quality measures such as estimation variance or mean squared error, 19,59 but in practice, these methods tend 
to be rather computationally demanding. Instead we simply multiply the data d(r ) by the Slepian functions or 
"tapers" g a (f ) designed for the region of interest prior to computing power and then averaging: 



(i+1)2 '4*- * ' 



ct—1 v / m— — l 



(71) 



Bias of the multitaper estimate 

The expected value of the estimate (71) is 

l+L 

(^ MT >= £ M W {S V +N V ), (72) 

l'=l-L 

where the eigenvalue- weighted multitaper coupling matrix, using Wigner 3-j functions, 60,61 is given by 

It is remarkable that this result depends only upon the chosen bandwidth L and is completely independent of 
the size, shape or connectivity of the region R, even as R = f2. Moreover, every row of the matrix in eq. (73) 
sums to unity, which ensures that a the multitaper spectral estimate S , i MT has no leakage bias in the case of a 
perfectly white spectrum provided the noise bias is subtracted as well: (5 , ; MT ) — MaiNf = S if Si = S. 

Variance of the multitaper estimate 

Under the moderately colored approximation, which is more easily justified in this case because the coupling (73) 
is confined to a narrow band of width less than or equal to 2L + 1, with L the bandwidth of the tapers, the 
eigenvalue- weighted multitaper covariance is 

^ = ^(S l +N l ){S v +N l ,)^(2p+l)T p ( l I [\\ (74) 

p=0 ^ ' 



where, using Wigner 3-j and 6-j functions, 60,61 
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In this expression B e , the boxcar power (64), which we note does depend on the shape of the region of interest, 
appears again, summed over angular degrees limited by 3-j selection rules to < e < 2L. The sum in eq. (74) 
is likewise limited to degrees < p < 2L. The effect of tapering with windows bandlimited to L is to introduce 
covariance between the estimates at any two different degrees I and I' that are separated by fewer than 2L + 1 
degrees. Eqs (74)-(75) are very efficiently computable, which should make them competitive with, e.g., jackknifed 
estimates of the estimation variance. 58,62,63 

The crux of the analysis lies in the fact that the matrix of the spectral covariances between sing/e-tapered 
estimates is almost diagonal, 19 showing that the individual estimates that enter the weighted formula (71) are 
almost uncorrelated statistically. This embodies the very essence of the multitaper method. It dramatically 
reduces the estimation variance at the cost of small increases of readily quantifiable bias. 

4. PRACTICAL CONSIDERATIONS 

In this section we now turn to the very practical context of sampled, e.g. geodetic, data on the sphere. We shall 
deal exclusively with bandlimited functions, which are equally well expressed in the spherical harmonic as the 
Slepian basis, namely: 

L l (L + lf 

/( f ) = E E fimYimir) - Yl f°9 a (P), (76) 

Z=0 m= — l a=l 

whereby the Slepian-basis expansion coefficients are obtained as 

f a = f f(v)g a (v)dn. (77) 

If the function of interest is spatially localized in the region R, a truncated reconstruction using Slepian functions 
built for the same region will constitute a very good, and sparse, local approximation to it: 64 

N 3D 

a=l 

We represent any sampled, bandlimited function / by an M-dimensional column vector 

f=(/i ••• fj ■■■ /m) T , (79) 

where fj = /(fj) is the value of / at pixel j, and M is the total number of sampling locations. In the most general 
case the distribution of pixel centers will be completely arbitrary. The special case of equal-area pixelization of a 
2-D function /(f) on the unit sphere is analogous to the equispaced digitization of a 1-D time series. Integrals 
will then be assumed to be approximated with sufficient accuracy by a Riemann sum over a dense set of pixels, 

/"/(f)dfiw AQ^fj and J f(r ) dtl w Afl f T f . (80) 



We have deliberately left the integration domain out of the above equations to cover both the cases of sampling 
over the entire unit sphere surface O, in which case the solid angle Ail = A-k/M (case 1) as well as over an 



incomplete subdomain R C 0, in which case AO = A/M, with A the area of the region R (case 2). If we collect 
the real spherical harmonic basis functions Yi m into an (L + l) 2 x M-dimcnsional matrix 



/ Foo(fi) ••• r o(f, 



ioo (?m) \ 



(81) 



and the spherical harmonic coefficients of the function into an (L + l) 2 x 1-dimcnsional vector 

f = ( /oo • • • /zm • • • III ) T , (82) 

we can write the spherical harmonic synthesis in eq. (76) for sampled data without loss of generality as 

f = Y T f . (83) 

We will adhere to the notation convention of using sans-serif fonts (e.g. f , Y) for vectors or matrices that depend 
on at least one spatial variable, and serifed fonts (e.g. f , D ) for those that are entirely composed of "spectral" 
quantities. In the case of dense, equal-area, whole-sphere sampling we have an approximation to eq. (29): 



YY T w AO" 1 ! (case 1), 



(84) 



where the elements of the (L + l) 2 x (L + l) 2 -dimcnsional spectral identity matrix I are given by the Kronecker 
deltas 5u'5 mm '. In the case of dense, equal-area, sampling over some closed region R, we find instead an 
approximation to the (L + l) 2 x (L + l) 2 -dimensional "spatiospectral localization matrix": 

YY T w AO _1 D (case 2), (85) 

where the elements of D arc those defined in eq. (33b). 

Let us now introduce the (L + l) 2 x (L + l) 2 -dimensional matrix of spectral Slepian eigenfunctions by 

/ ffooi ••• flooa ••• 5oo(L+i) 2 \ 



G 



9lr. 



(86) 



\ 9LL1 ■■■ 9LL a ■■■ 9LL{L+l) 2 J 

This is the matrix that contains the eigenfunctions of the problem defined in eq. (33), which we rewrite as 

DG = G A, (87) 

where the diagonal matrix with the concentration eigenvalues is given by 

A = diag ( Ai • • • \ a ■■■ A (L+1)2 ) . (88) 

The spectral orthogonality relations of eq. (35) are 

G T G = I, G T DG = A, (89) 

where the elements of the (L + l) 2 x (L + l) 2 -dimensional Slepian identity matrix I are given by the Kronecker 
deltas 5 a p- We write the Slepian functions of eq. (31) as 



G = G T Y and Y = GG, 



(90) 



where the (L + l) 2 x M-dimensional matrix holding the sampled spatial Slepian functions is given by 



G ••• g a (ij) ■■■ ■ (91) 

V 3(l+i)2(£i) ••• 5(L+i) 2 (fj) ••• 9(L+iy{r M ) ) 

Under a dense, equal-area, whole-sphere sampling, we will recover the the spatial orthogonality of cq. (36) 
approximately as 

GG T wAfi- 1 I (easel), (92) 
whereas for dense, equal-area, sampling over a region R we will get, instead, 

GG^Afl^A (case 2). (93) 

With this matrix notation we shall revisit both estimation problems of the previous section. 

4.1. Problem (i), revisited 
Spherical harmonic solution 

If we treat eq. (83) as a noiseless inverse problem in which the sampled data f are given but from which the 
coefficients f arc to be determined, we find that for dense, equal-area, whole-sphere sampling, the solution 

fwAOYf (easel) (94) 

is simply the discrete approximation to the spherical harmonic analysis formula (25). For dense, equal-area, 
regional sampling we need to calculate 

fsiASJD-'Yf (case 2). (95) 

Both of these cases are simply the relevant solutions to the familiar overdetermined spherical harmonic inversion 
problem 65-67 for discretely sampled data, i.e. the least-squares solution to eq. (83), 

f=(YY T )- 1 Yf, (96) 

for the particular cases described by eqs (84)-(85). In eq. (95) we furthermore recognize the discrete version 
of eq. (57) with r\ = 0, the undamped solution to the minimum mean squared error inverse problem posed in 
continuous form in eq. (56). From the continuous limiting case eq. (57) we thus discover the general form that 
damping should take in regularizing the ill-conditioned inverse required in eqs (95)-(96). Its principal property 
is that it differs from the customary ad hoc practice of adding small values on the diagonal only. Finally, in 
the most general and admittedly most commonly encountered case of randomly scattered data we require the 
Moore-Penrose pseudo-inverse 

f = pinv(Y T )f, (97) 

which is constructed by inverting the singular value decomposition (svd) of Y T with its singular values truncated 
beyond where they become vanishingly small. 68 Solving eq. (97) by truncated svd is equivalent to inverting a 
truncated eigenvalue expansion of the normal matrix YY T as it appears in eq. (96), as can be easily shown. 

Slepian basis solution 

If we collect the Slepian expansion coefficients of the function / into the (L + l) 2 x 1-dimensional vector 

t = (A ••• f a ••• f( L+ ir) T , (98) 

the expansion (76) in the Slepian basis takes the form 

f = G T t = Y T Gt, (99) 



where we used eqs (89)-(90) to obtain the second equality. Comparing eq. (99) with eq. (83), we see that the 
Slepian expansion coefficients of a function transform to and from the spherical harmonic coefficients as: 



f = Gt and t = G T f. (100) 

Under dense, equal-area, sampling with complete coverage the coefficients in eq. (99) can be estimated from 

UAfiGf (easel), (101) 

the discrete, approximate, version of eq. (77). For dense, equal-area, sampling in a limited region R we get 

t w AO A _1 Gf (case 2). (102) 

As expected, both of the solutions (101)-(102) are again special cases of the overdetermined least-squares solution 

t = (GG T ) _1 Gf, (103) 

as applied to eqs (92)-(93). We encountered eq. (102) before in the continuous form of eq. (59); it solves 
the undamped minimum mean squared error problem (56). The regularization of this ill-conditioned inverse 
problem may be achieved by truncation of the concentration eigenvalues, e.g. by restricting the size of the 
(L + l) 2 x (L + l) 2 -dimensional operator GG T to its first 7V 3D x A 3D subblock. Finally, in the most general, 
scattered-data case, we would be using an eigenvalue-truncated version of eq. (103), or, which is equivalent, form 
the pseudo-inverse 

t = pinv(G T )f. (104) 

The solutions (94)-(96) and (101)— (103) are equivalent and differ only by the orthonormal change of basis 
from the spherical harmonics to the Slepian functions. Indeed, using eqs (100) and (90) to transform eq. (103) 
into an equation for the spherical harmonic coefficients, and comparing with eq. (96) exposes the relation 

G(GG T ) _1 G T = (YY T )-\ (105) 

which is a trivial identity for case 1 (insert eqs 84, 92 and 89) and, after substituting eqs (85) and (93), entails 

GA _1 G T = D 1 (106) 

for case 2, which holds by virtue of eq. (89). Eq. (105) can also be verified directly from eq. (90), which implies 

YY T = G(GG T )G T . (107) 

The popular but labor-intensive procedure by which the unknown spherical harmonic expansion coefficients of a 
scattered data set are obtained by forming the Moore-Penrose pseudo-inverse as in eq. (97) is thus equivalent to 
determining the truncated Slepian solution of eq. (102) in the limit of continuous and equal-area, but incomplete 
data coverage. In that limit, the generic eigenvalue decomposition of the normal matrix becomes a specific 
statement of the Slepian problem as we encountered it before, namely 

YY T AO = US 2 U T -» D = GAG T . (108) 

Such a connection has been previously pointed out for time series 69 and leads to the notion of "generalized prolate 
spheroidal functions" 70 should the "Slepian" functions be computed from a formulation of the concentration 
problem in the scattered data space directly, rather than being determined by sampling those obtained from 
solving the corresponding continuous problem, as we have described here. 

Above, we showed how to stabilize the inverse problem of eq. (96) by damping. We dealt with the case of 
continuously available data only; the form in which it appears in eq. (57) makes it clear that damping is hardly 
practical for scattered data. Indeed, it requires knowledge of the complementary localization operator D, in 
addition to being sensitive to the choice of r], whose optimal value depends implicitly on the unknown signal-to- 
noise ratio. 17 The data-driven approach taken in eq. (97) is the more sensible one. 68 We have now seen that, 



in the limit of continuous partial coverage, this corresponds to the optimal solution of the problem formulated 
directly in the Slepian basis. It is consequently advantageous to also work in the Slepian basis in case the data 
collected are scattered but closely collocated in some region of interest. Prior knowledge of the geometry of this 
region and a prior idea of the spherical harmonic bandwidth of the data to be inverted allows us to construct 
a Slepian basis for the situation at hand, and the problem of finding the Slepian expansion coefficients of the 
unknown underlying function can be solved using eqs (103)-(104). The measure within which this approach 
agrees with the theoretical form of eq. (102) will depend on how regularly the data are distributed within 
the region of study, i.e on the error in the approximation (93). But if indeed the scattered-data Slepian normal 
matrix GG T is nearly diagonal in its first iV 3D x iV 3D -dimensional block due to the collocated observations having 
been favorably, if irregularly, distributed, then eq. (102), which, strictly speaking, requires no matrix inversion, 
can be applied directly. If this is not the case, but the data are still collocated or we are only interested in a 
local approximation to the unknown signal, we can restrict G to its first N 3D rows, prior to diagonalizing GG T or 
performing the svd of a partial G T as necessary to calculate eqs (103)-(104). Compared to solving eqs (96)-(97), 
the computational savings will still be substantial, as only when R w Q, will the operator YY T be nearly diagonal. 
Truncation of the eigenvalues of YY T is akin to truncating the matrix GG T itself, which is diagonal or will be 
nearly so. With the theoretically, continuously determined, sampled Slepian functions as a parametrization, the 
truncated expansion is easy to obtain and the solution will be locally faithful within the region of interest R. 
In contrast, should we truncate YY T itself, without first diagonalizing it, we would be estimating a low-degree 
approximation of the signal which would have poor resolution everywhere. 

Bias and variance 

For completeness we briefly return to the expressions for the mean squared estimation error of the damped 
spherical-harmonic and the truncated Slepian function methods, eqs. (61)-(62), which we quoted for the example 
of "white" signal and noise. Introducing the (L + l) 2 x (L + l) 2 -dimcnsional spectral matrices 

H = A + t?(I-A), (109a) 

V = iVH- 2 A, and B = \/£H _1 (I — A), (109b) 

we handily rewrite the "full" version of eq. (61) in two spatial variables as the error covariance matrix 

(e(f)e(f')) =G T (V + 77 2 B 2 )G. (110) 

We subdivide the matrix with Slepian functions into the truncated set of the best-concentrated a = 1 — ► iV 3D 
and the complementary set of remaining a = N 3D + 1 — > (L + l) 2 functions, as follows 

G = ( G T G T ) T , (111) 

and similarly separate the eigenvalues, writing 

A = diag(Ai ••• X N 3n ), (112a) 
A = diag( Ajv3d +1 ••• A (i+1) 2) . (112b) 

Likewise, the identity matrix is split into two parts, I and I. If we now also redefine 

V^iVA -1 , and B = VSI, (113a) 

Y = A/A -1 , and B = VS1, (113b) 
the equivalent version of eq. (62) is readily transformed into the full spatial error covariance matrix 

(e(r)e(f )) = G T YG + G T B 2 G. (114) 

In selecting the Slepian basis we have thus successfully separated the effect of the variance and the bias on the 
mean squared reconstruction error of a noisily observed signal. If the region of observation is a contiguous closed 
domain R C and the truncation takes place at the Shannon number iV 3D , we have thereby identified the 
variance as due to noise in the region where data are available, and the bias to signal neglected in the truncated 
expansion — which, in the proper Slepian basis, corresponds to the regions over which no observations exist. 

Finally, we shall also apply the notions of discretely acquired data to the solutions of problem (ii), below. 



4.2. Problem (ii), revisited 

We need two more pieces of notation in order to rewrite the expressions for the spectral estimates (65) and (71) in 
the "pixel-basis" . First we construct the M x M-dimensional symmetric spatial matrix collecting the fixed-degree 
Legendre polynomials evaluated at the angular distances between all pairs of observations points, 
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(115) 
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The elements of P/ are thus J2m=-i Yi m {*i)Ylm{* j), by the addition theorem, eq. (30). And finally, we define G" 
the M x M symmetric matrix with elements given by 
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(116) 



The spherical periodogram 

The expression equivalent to eq. (65) is now written as 
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whereby the column vector d contains the sampled data as in the notation for eq. (79). This lends itself easily 
to computation, and the statistics of eqs (66)-(69) hold, approximately, for sufficiently densely sampled data. 

The spherical multitaper estimate 

Finally, the expression equivalent to eq. (71) becomes 
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(118) 



Both eqs (117) and (118) are quadratic forms, earning them the nickname "quadratic spectral estimators". 71 
The key difference with the maximum-likelihood estimator popular in cosmology, 45-47 which can also be written 
as a quadratic form, 43 is that neither P; nor G" depend on the unknown spectrum itself. In contrast, maximum- 
likelihood estimation is inherently non-linear, requiring iteration to converge to the most probable estimate of 
the power spectral density. 20 

The estimate (118) has the statistical properties that we listed earlier as eqs (72)— (75). These continue to hold 
when the data pixelization is fine enough to have integral expressions of the kind (80) be exact. As mentioned 
before, for completely irregularly and potentially non-densely distributed discrete data on the sphere, "gener- 
alized" Slepian functions 70 could be constructed specifically for the purpose of their power spectral estimation, 
and used to build the operator (116). 



5. CONCLUSIONS 

What is the information contained in a bandlimited set of scientific observations made over an incomplete, e.g. 
temporally or spatially limited sampling domain? How can this "information", e.g. an estimate of the signal 
itself, or of its energy density, be determined from noisy data, and how shall it be represented? These seemingly 
age-old fundamental questions, which have implications beyond the scientific, 72 had been solved — some say, 
by conveniently ignoring them — heuristically, by engineers, well before receiving their first satisfactory answers 
given in the theoretical treatment by Slepian, Landau and Pollak; 2-4 first for "continuous" time series, later 



generalized to the multidimensional and discrete cases. 23 ' 32 ' By the "Slepian functions" in the title of this 
contribution, we have lumped together all functions that are "spatiospcctrally" concentrated, quadratically, in 
the original sense of Slepian. In one dimension, these are the "prolate spheroidal functions" whose popularity is 
as enduring as their utility. In two Cartesian dimensions, and on the surface of the unit sphere, their time is yet 
to come, and we have made a case for it in this study. 

The answers to the questions posed above are as ever relevant for the geosciences of today. There, we often face 
the additional complications of irregularly shaped study domains, scattered observations of noise-contaminated 
potential fields, perhaps collected from an altitude above the source by airplanes or satellites, and an acquisition 
and model-space geometry that is rarely if ever non-spherical. Thus the Slepian functions are especially suited 
for geoscientific applications and to study any type of geographical information, in general. 

Two problems that are of particular interest in the geosciences, but also further afield, are how to form a 
statistically "optimal" estimate of the signal giving rise to the data, and how to estimate the power spectral 
density of such signal. The first, an inverse problem that is linear in the data, applies to forming mass flux 
estimates from time-variable gravity, e.g. by the GRACE mission, or to the characterization of the terrestrial 
magnetic field by satellites such as CHAMP or swarm. The second, which is quadratic in the data, is of interest 
in studying the statistics of the Earth's or planetary topography, and especially for the cross-spectral analysis of 
gravity and topography, which can yield important clues about the internal structure of the planets. The second 
problem is also of great interest in cosmology, where missions such as wmap and PLANCK are mapping the cosmic 
microwave background radiation, which is best modeled spectrally to constrain models of the evolution of our 
universe. 

Slepian functions, as we have shown by focusing on the case of spherical geometry, provide the mathematical 
framework to solve such problems. They are a convenient and easily obtained doubly-orthogonal mathematical 
basis in which to express, and thus by which to recover, signals that are geographically localized, or incompletely 
(and noisily) observed. For this they are much better suited than the traditional Fourier or spherical harmonic 
bases, and they are more "geologically intuitive" than wavelet bases in retaining a firm geographic footprint and 
preserving the traditional notions of frequency or spherical harmonic degree. They are furthermore extremely 
performant as data tapers to regularize the inverse problem of power spectral density determination from noisy 
and patchy observations, which can then be solved satisfactorily without costly iteration. Finally, by the inter- 
pretation of the Slepian functions as their limiting cases, much can be learned about the statistical nature of 
such inverse problems when the data provided are themselves scattered within a specific areal region of study. 
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